home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / WINDOWS / WXLSLIB.ARJ / STATS.LSP < prev    next >
Lisp/Scheme  |  1992-02-20  |  11KB  |  325 lines
  1. ;;;;
  2. ;;;; statistics.lsp XLISP-STAT statistics functions
  3. ;;;; XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney
  4. ;;;; Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz
  5. ;;;; You may give out copies of this software; for conditions see the file
  6. ;;;; COPYING included with this distribution.
  7. ;;;;
  8.  
  9. (provide "statistics")
  10.  
  11. ;;;;
  12. ;;;; Data File Reading 
  13. ;;;;
  14.  
  15. (defun count-file-columns (fname)
  16. "Args: (fname)
  17. Returns the number of lisp items on the first nonblank line of file FNAME."
  18.   (with-open-file (f fname)
  19.     (if f
  20.         (let ((line (do ((line (read-line f) (read-line f))) 
  21.                         ((or (null line) (< 0 (length line))) line))))
  22.           (if line
  23.               (with-input-from-string (s line)
  24.                 (do ((n 0 (+ n 1)) (eof (gensym))) 
  25.                     ((eq eof (read s nil eof)) n))))))))
  26.  
  27. (defvar *xlisptable* *readtable*)
  28.  
  29. (if (not (fboundp 'open-file-dialog))
  30.   #+dialogs
  31.   (defun open-file-dialog (&optional set)
  32.     (get-string-dialog "Enter a data file name:"))
  33.   #-dialogs
  34.   (defun open-file-dialog (&optional set)
  35.     (error "You must provide a file name explicitly")))
  36.  
  37. (defun read-data-file (&optional (file (open-file-dialog t)))
  38. "Args:  (file)
  39. Returns a list of all lisp objects in FILE. FILE can be a string or a symbol,
  40. in which case the symbol'f print name is used."
  41.   (if file
  42.       (let ((oldtable *readtable*)
  43.             (oldbreak *breakenable*)
  44.             (eof (gensym)))
  45.         (setq *readtable* *xlisptable*)
  46.         (setq *breakenable* nil)
  47.         (with-open-file (f file)
  48.           (if f
  49.               (unwind-protect
  50.                (do* ((r (read f nil eof) (read f nil eof))
  51.                      (x (list nil))
  52.                      (tail x (cdr tail)))
  53.                     ((eq r eof) (cdr x))
  54.                     (setf (cdr tail) (list r)))
  55.                (setq *breakenable* oldbreak)
  56.                (setq *readtable* oldtable)))))))
  57.  
  58. ;;; New definition to avoid stack size limit in apply
  59. (defun read-data-columns (&optional (file (open-file-dialog t))
  60.                                     (cols (if file 
  61.                                               (count-file-columns file))))
  62. "Args: (&optional file cols)
  63. Reads the data in FILE as COLS columns and returns a list of lists representing the columns."
  64.   (if (and file cols)
  65.       (transpose (split-list (read-data-file file) cols))))
  66.  
  67. #+unix
  68. (defun load-data (file)
  69. "Args: (file)
  70. Read in data file from the data examples library."
  71.   (if (load (format nil "~aData/~a" *default-path* file))
  72.       t
  73.       (load (format nil "~aExamples/~a" *default-path* file))))
  74.  
  75. #+unix
  76. (defun load-example (file)
  77. "Args: (file)
  78. Read in lisp example file from the examples library."
  79.   (if (load (format nil "~aExamples/~a" *default-path* file))
  80.       t
  81.       (load (format nil "~aData/~a" *default-path* file))))
  82. #+macintosh
  83. (defun load-data (s) (require s (concatenate 'string ":Data:" s)))
  84. #+macintosh
  85. (defun load-example (s) (require s (concatenate 'string ":Examples:" s)))
  86.  
  87. #+msdos
  88. (defun load-data (file)
  89. "Args: (file)
  90. Read in data file from the data examples library."
  91.   (load (format nil "~aData\\~a" *default-path* file)))
  92.  
  93. #+msdos
  94. (defun load-example (file)
  95. "Args: (file)
  96. Read in lisp example file from the examples library."
  97.   (load (format nil "~aExamples\\~a" *default-path* file)))
  98.  
  99. ;;;;
  100. ;;;; Listing and Saving Variables and Functions
  101. ;;;;
  102.  
  103. (defvar *variables*)
  104. (defvar *ask-on-redefine*)
  105.  
  106. (defmacro def (symbol value)
  107. "Syntax: (def var form)
  108. VAR is not evaluated and must be a symbol.  Assigns the value of FORM to
  109. VAR and adds VAR to the list *VARIABLES* of def'ed variables. Returns VAR.
  110. If VAR is already bound and the global variable *ASK-ON-REDEFINE*
  111. is not nil then you are asked if you want to redefine the variable."
  112.   `(unless (and *ask-on-redefine*
  113.                 (boundp ',symbol)
  114.                 (not (y-or-n-p "Variable has a value. Redefine?")))
  115.            (pushnew ',symbol *variables*)
  116.            (setf ,symbol ,value)
  117.            ',symbol))
  118.   
  119. (defun variables-list () 
  120.     (mapcar #'intern (sort-data (mapcar #'string *variables*))))
  121.  
  122. (defun variables ()
  123. "Args:()
  124. Returns a list of the names of all def'ed variables to STREAM"
  125.   (if *variables*
  126.       (mapcar #'intern (sort-data (mapcar #'string *variables*)))))
  127.   
  128. (defun savevar (vars file)
  129. "Args: (vars file-name-root)
  130. VARS is a symbol or a list of symbols. FILE-NAME-ROOT is a string (or a symbol
  131. whose print name is used) not endinf in .lsp. The VARS and their current values
  132. are written to the file FILE-NAME-ROOT.lsp in a form suitable for use with the
  133. load command."
  134.   (let ((f (open (strcat (string file) ".lsp") :direction :output))
  135.         (vars (if (consp vars) vars (list vars)))
  136.         (oldbreak *breakenable*))
  137.     (setq *breakenable* nil)
  138.     (unwind-protect
  139.       (mapcar
  140.         (lambda (x)
  141.             (if (objectp (symbol-value x))
  142.                 (print `(def ,x ,(send (symbol-value x) :save)) f)
  143.                 (print `(def ,x ',(symbol-value x)) f)))
  144.         vars)
  145.       (setq *breakenable* oldbreak)
  146.       (close f))
  147.     vars))
  148.  
  149. (defun undef (v)
  150. "Args: (v)
  151. If V is the symbol of a defined variable the variable it is unbound and
  152. removed from the list of defined variables. If V is a list of variable
  153. names each is unbound and removed. Returns V."
  154.   (dolist (s (if (listp v) v (list v)))
  155.           (when (member s *variables*)
  156.                 (setq *variables* (delete s *variables*))
  157.                 (makunbound s)))
  158.   v)
  159.         
  160.  
  161. ;;;;
  162. ;;;; Basic Summary Statistics
  163. ;;;;
  164.  
  165. (defun standard-deviation (x)
  166. "Args: (x)
  167. Returns the standard deviation of the elements x. Vector reducing."
  168.   (let ((n (count-elements x))
  169.         (r (- x (mean x))))
  170.     (sqrt (* (mean (* r r)) (/ n (- n 1))))))
  171.  
  172. (defun quantile (x p)
  173. "Args: (x p)
  174. Returns the P-th quantile(s) of sequence X. P can be a number or a sequence."
  175.   (let* ((x (sort-data x))
  176.          (n (length x))
  177.          (np (* p (- n 1)))
  178.          (low (floor np))
  179.          (high (ceiling np)))
  180.     (/ (+ (select x low) (select x high)) 2)))
  181.     
  182. (defun median (x) 
  183. "Args: (x)
  184. Returns the median of the elements of X."
  185.   (quantile x 0.5))
  186.  
  187. (defun interquartile-range (x) 
  188. "Args: (number-data)
  189. Returns the interquartile range of the elements of X."
  190.   (apply #'- (quantile x '(0.75 0.25))))
  191.  
  192. (defun fivnum (x) 
  193. "Args: (number-data)
  194. Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
  195. max) of the elements X."
  196.   (quantile x '(0 .25 .5 .75 1)))
  197.  
  198. (defun covariance-matrix (&rest args)
  199. "Args: (&rest args)
  200. Returns the sample covariance matrix of the data columns in ARGS. ARGS may
  201. consist of lists, vectors or matrices."
  202.   (let ((columns (apply #'append 
  203.                         (mapcar (lambda (x) 
  204.                                   (if (matrixp x) (column-list x) (list x)))
  205.                                 args))))
  206.     (/ (cross-product (apply #'bind-columns 
  207.                              (- columns (mapcar #'mean columns))))
  208.        (- (length (car columns)) 1))))
  209.  
  210. ;;;;
  211. ;;;; Basic Sequence Operations
  212. ;;;;
  213.  
  214. (defun difference (x)
  215. "Args: (x)
  216. Returns differences for a sequence X."
  217.   (let ((n (length x)))
  218.     (- (select x (iseq 1 (1- n))) (select x (iseq 0 (- n 2))))))
  219.  
  220. (defun rseq (a b num)
  221. "Args: (a b num)
  222. Returns a list of NUM equally spaced points starting at A and ending at B."
  223.   (+ a (* (iseq 0 (1- num)) (/ (- b a) (1- num)))))
  224.  
  225.  
  226. ;;;;
  227. ;;;; Linear Algebra Functions
  228. ;;;;
  229.  
  230. (defun matrix (dim data)
  231. "Args: (dim data)
  232. returns a matrix of dimensions DIM initialized using sequence DATA
  233. in row major order." 
  234.   (let ((dim (coerce dim 'list))
  235.         (data (coerce data 'list)))
  236.     (make-array dim :initial-contents (split-list data (nth 1 dim)))))
  237.  
  238. (defun print-matrix (a &optional (stream *standard-output*))
  239. "Args: (matrix &optional stream)
  240. Prints MATRIX to STREAM in a nice form that is still machine readable"
  241.   (unless (matrixp a) (error "not a matrix - ~a" a))
  242.   (let ((size (min 15 (max (map-elements #'flatsize a)))))
  243.     (format stream "#2a(~%")
  244.     (dolist (x (row-list a))
  245.             (format stream "    (")
  246.             (let ((n (length x)))
  247.               (dotimes (i n)
  248.                        (let ((y (aref x i)))
  249.                          (cond
  250.                            ((integerp y) (format stream "~vd" size y))
  251.                            ((floatp y) (format stream "~vg" size y))
  252.                            (t (format stream "~va" size y))))
  253.                        (if (< i (- n 1)) (format stream " "))))
  254.             (format stream ")~%"))
  255.     (format stream "   )~%")
  256.     nil))
  257.  
  258. (defun solve (a b)
  259. "Args: (a b)
  260. Solves A x = B using LU decomposition and backsolving. B can be a sequence
  261. or a matrix."
  262.   (let ((lu (lu-decomp a)))
  263.     (if (matrixp b)
  264.         (apply #'bind-columns 
  265.                (mapcar #'(lambda (x) (lu-solve lu x)) (column-list b)))
  266.         (lu-solve lu b))))
  267.         
  268. (defun backsolve (a b)
  269. "Args: (a b)
  270. Solves A x = B by backsolving, assuming A is upper triangular. B must be a
  271. sequence. For use with qr-decomp."
  272.   (let* ((n (length b))
  273.          (sol (make-array n)))
  274.     (dotimes (i n)
  275.              (let* ((k (- n i 1))
  276.                     (val (elt b k)))
  277.                (dotimes (j i)
  278.                         (let ((l (- n j 1)))
  279.                           (setq val (- val (* (aref sol l) (aref a k l))))))
  280.                (setf (aref sol k) (/ val (aref a k k)))))
  281.     (if (listp b) (coerce sol 'list) sol)))
  282.  
  283. (defun eigenvalues (a) 
  284. "Args: (a)
  285. Returns list of eigenvalues of square, symmetric matrix A"
  286.   (first (eigen a)))
  287.  
  288. (defun eigenvectors (a) 
  289. "Args: (a)
  290. Returns list of eigenvectors of square, symmetric matrix A"
  291.   (second (eigen a)))
  292.  
  293. (defun accumulate (f s)
  294. "Args: (f s)
  295. Accumulates elements of sequence S using binary function F.
  296. (accumulate #'+ x) returns the cumulative sum of x."
  297.   (let* ((result (list (elt s 0)))
  298.          (tail result))
  299.     (flet ((acc (dummy x)
  300.                 (rplacd tail (list (funcall f (first tail) x)))
  301.                 (setf tail (cdr tail))))
  302.       (reduce #'acc s))
  303.     (if (vectorp s) (coerce result 'vector) result)))
  304.  
  305. (defun cumsum (x)
  306. "Args: (x)
  307. Returns the cumulative sum of X."
  308.   (accumulate #'+ x))
  309.  
  310. (defun combine (&rest args) 
  311. "Args (&rest args) 
  312. Returns sequence of elements of all arguments."
  313.   (copy-seq (element-seq args)))
  314.  
  315. (defun lowess (x y &key (f .25) (steps 2) (delta -1) sorted)
  316. "Args: (x y &key (f .25) (steps 2) delta sorted)
  317. Returns (list X YS) with YS the LOWESS fit. F is the fraction of data used for
  318. each point, STEPS is the number of robust iterations. Fits for points within
  319. DELTA of each other are interpolated linearly. If the X values setting SORTED
  320. to T speeds up the computation."
  321.   (let ((x (if sorted x (sort-data x)))
  322.         (y (if sorted y (select y (order x))))
  323.         (delta (if (> delta 0.0) delta (/ (- (max x) (min x)) 50))))
  324.     (list x (|base-lowess| x y f steps delta))))
  325.